clear all;
clf;
clc;
close all;

fs=10;NN=1270;    %采样频率和数据点数；采样定理，频域最大的频率绝对值是时域采样间隔倒数的一半

n=0:NN-1;t=n/fs;  %时间序列
f_p=1;
f_c=10e6;
%----PRBS---Start--------------------------------
N = 127; % lenght of the PRBS sequence
M = 10; % number of samples per bit
uk = idinput([N 1 100],'prbs',[0 1],[-1 1]);  %u = idinput([Period,Nu,NumPeriod]) 返回具有指定周期和周期数的 Nu 通道周期性随机二进制输入信号。 每个输入通道信号的长度为 NumPeriod*Period。
if M<1 
    fprintf('\t\t ERROR\t\t\n');
    fprintf('Number of samples per bit is less than 1\n');
    return;
end
for k=1:N
    if (uk(k) == -1)
        for i=1:M
        waveformTemp(k,i) = [-1];
        end
        waveform{k} = waveformTemp(k,:);
    elseif (uk(k) == 1)
        for i=1:M
        waveformTemp(k,i) = [1];
        end
        waveform{k} = waveformTemp(k,:);
    end  
end
PRBSwave = cell2mat(waveform);
figure;
h = plot(PRBSwave,'k');
set (h, 'LineWidth', 2); % for a width of n
axis([0,length(PRBSwave),-1.5,1.5]);
title('PRBS waveform');
xlabel('samples');
ylabel('PRBS +1 / -1 values');
%%----PRBS---End---------------------------------
%-----------PRBS产生方法二------------------%
% n=7; % 阶次 
% p=2^n-1; % 循环周期 
% pr = idinput([p 1 20], 'prbs',[0,1]);
% PRBSwave=pr';
% stairs(PRBSwave);
%---频率f的定义-2-开始-----%
pm=0.5*PRBSwave;
beta=1.*pi./2;

f=(-NN/2:NN/2-1)*1*fs/NN;

y=exp(1i.*(2.*pi.*f_c.*t+beta.*pm))+randn(size(t))/10; %信号
y1=exp(1i.*(4.*pi.*f_c.*t+2*beta.*pm))+randn(size(t))/10; %信号

y2=exp(1i.*(2.*pi.*f_c.*t+2*beta.*pm))+randn(size(t))/10; %信号
y3=exp(1i.*(4.*pi.*f_c.*t+4*beta.*pm))+randn(size(t))/10; %信号

y4=exp(1i.*(2.*pi.*f_c.*t))+randn(size(t))/10; %信号
y5=exp(1i.*(4.*pi.*f_c.*t))+randn(size(t))/10; %信号

S_PM=fftshift(fft(y,length(t)))./length(t);
S_PM1=fftshift(fft(y1,length(t)))./length(t);

S_PM2=fftshift(fft(y2,length(t)))./length(t);
S_PM3=fftshift(fft(y3,length(t)))./length(t);

S_PM4=fftshift(fft(y4,length(t)))./length(t);
S_PM5=fftshift(fft(y5,length(t)))./length(t);

I_PM=abs(S_PM).^2;
I_PM1=abs(S_PM1).^2;
I_PM2=abs(S_PM2).^2;
I_PM3=abs(S_PM3).^2;
I_PM4=abs(S_PM4).^2;
I_PM5=abs(S_PM5).^2;

figure;
plot(f,I_PM);   %对信号进行快速Fourier变换,求得Fourier变换后的振幅
% plot(f,20*log10(abs(S_PM)));
title(['Spectrum of Fundamental Laser {\gamma}=',num2str(beta)],'FontName','Times New Roman','FontWeight','Bold');
xlabel('Normalized Frequency');ylabel('Spectral Power [a. u.]');

figure;
plot(f,I_PM1,'color','k');   %对信号进行快速Fourier变换,求得Fourier变换后的振幅
% plot(f,20*log10(abs(S_PM1)));
title(['Spectrum of Second Harmonic {\gamma}={\pi}'],'FontName','Times New Roman','FontWeight','Bold');
xlabel('Normalized Frequency');ylabel('Spectral Power [a. u.]');

figure;
plot(f,I_PM2);   %对信号进行快速Fourier变换,求得Fourier变换后的振幅
% plot(f,20*log10(abs(S_PM2)));
title(['Spectrum of Fundamental Laser {\gamma}={\pi}'],'FontName','Times New Roman','FontWeight','Bold');
xlabel('Normalized Frequency');ylabel('Spectral Power [a. u.]');

figure;
plot(f,I_PM3,'color','k');   %对信号进行快速Fourier变换,求得Fourier变换后的振幅
% plot(f,20*log10(abs(S_PM3)));
title(['Spectrum of Second Harmonic {\gamma}=2{\pi}'],'FontName','Times New Roman','FontWeight','Bold');
xlabel('Normalized Frequency');ylabel('Spectral Power [a. u.]');

figure;
plot(f,I_PM4);   %对信号进行快速Fourier变换,求得Fourier变换后的振幅
% plot(f,20*log10(abs(S_PM4)));
title(['Spectrum of Fundamental Laser {\gamma}=0'],'FontName','Times New Roman','FontWeight','Bold');
xlabel('Normalized Frequency');ylabel('Spectral Power [a. u.]');

figure;
plot(f,I_PM5,'color','k');   %对信号进行快速Fourier变换,求得Fourier变换后的振幅
% plot(f,20*log10(abs(S_PM5)));
title(['Spectrum of Second Harmonic {\gamma}=0'],'FontName','Times New Roman','FontWeight','Bold');
xlabel('Normalized Frequency');ylabel('Spectral Power [a. u.]');


